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Abstract. We study the evolution with redshift of the num- 
ber density of Lya clouds in a CDM model using numeri- 
cal simulations including photo-ionization and cooling of the 
baryonic component. The photo-ionizing flux is consistently 
taken as proportional to the rate at which material cools below 
5000 K in the simulation. The number density of lines with 
l/?j(HI) > 10 14 cm~ 2 given by the simulation can be approxi- 
mated as AnIAz ~ 16.8(l+^)°- 8 +0.687(l+z) 2 - 5 which fits well 
^JjjE)data at any redshift. The number of weaker lines is repro- 
flBced within a factor of two and does not decrease with time 
•^fast as the number of strong lines. For A^(HI) > 10 13 cm~ 2 , 
find AnIAz oc (1+z) - 7 . 

At z = 0, the two-point correlation function for Lymana 
clouds with iV(HI) > 10 13 and > 10 14 cm~ 2 has a slope ~ - 
^^5 comparable to the slope of the two point correlation func- 
(tiqn for the galaxy distribution. The signal is higher for stronger 
rfetes for separations smaller than 1 Mpc. Part of the Lymana 
QQ-est is indeed associated with galaxies however the overall 
tture is that Lya absorption lines originate in the warm gas 
at traces the potential wells of filamentary structures defined 
: the dark matter. 

i^ey words: large-scale structure, intergalactic medium, 
^asars: absorption lines 
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1. Introduction 

Observation of QSO absorption line systems is the most sen- 
sitive way to probe baryonic material in the whole redshift 
range < z < 5. Among absorption line systems, Lya sys- 
tems without detected associated metal lines are the most nu- 
merous (about 150 per unit redshift at z ~ 2). They have been 
extensively studied, since they may trace the evolution of the in- 
tergalactic medium and the UV background. Indeed the clouds 
could be an intergalactic population possibly of primordial gas 
evolving independently of galaxies under the direct influence 
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of the intergalactic medium, the ionizing background, and the 
gravitational potential of dark matter. They are observed over 
a large redshift range, and they thus probe the evolution of the 
baryonic material and the structures formed by dark matter. 

The correct picture for the origin of the Lya forest will most 
likely involve a mixture of all the ideas invoked up to now, and 
probably the Lya forest arises through clouds in very different 
environments. Gas weakly enriched in heavy elements could be 
correlated with galaxies due to a variety of mechanisms: First, 
when galaxies form, the energy released could eject a substantial 
amount of gas at high velocity and thus to a comparatively large 
distance from the center (Lake 1988). This gas could be weakly 
enriched in heavy elements. Second, after formation, an initial 
starburst phase would inevitably lead to supernovae, creating a 
strong temperature gradient, ejecting metal-enriched gas, and 
thus again enhancing the medium surrounding the galaxy with 
metals. While large galaxies would be expected to undergo nu- 
merous starburst phases, resulting in roughly solar abundances, 
dwarf galaxies might be expected to have only one significant 
starburst phase, resulting in a weak metal enhancement in sur- 
rounding clouds. Third, during the process of merging, which 
seems to be a common event in the course of galaxy evolution, 
metal enriched gas could be ejected very far from the center 
(Morris & van den Bergh 1994), but still close enough to be 
associated with galaxies and thus show some clustering. 

However, primordial gas should be anticorrelated with star 
forming regions. The subpopulation of intergalactic clouds with 
primordial abundances should evolve independently as pres- 
sure confined gas (Ikeuchi & Ostriker 1986), or gravitation- 
ally confined gas in the potential wells of mini-halos (Ikeuchi 
1986, Rees 1986), such as those predicted in a cold dark matter 
(CDM) scenario. The model of gravitational collapse on sub- 
galactic scales (Bond et al. 1988) has received some support 
from recent hydrodynamic simulations (Cen et al. 1994) sug- 
gesting that Lya clouds could result from pancake structures 
(Zel'dovich 1970) in which the photoionized gas is maintained 
by a balance between gravitation and ram pressure due to in- 
falling gas. 

It is therefore important to clarify the link between absorp- 
tion line systems and galaxies. This can be done by studying 



the content of the gas in metals, the clustering properties of the 
systems, and the properties of the galaxies close to the line of 
sight to the quasar. 

It is generally assumed that the metallicity is very small 
(Chaffe et al. 1986), but this assumption is not firmly estab- 
lished. In fact, the presence of CIV in Lya systems has been 
claimed by Lu (1991), who finds evidence for the presence 
of the OVA 1548 line in a composite spectrum obtained by 
shifting and stacking 172 Lya systems of equivalent widths 
0.4 <W Y < 0.8 X. A similar study but using higher resolution 
data leads to an upper limit of one-hundredth of solar abun- 
dances (Tytler & Fan 1994). Very recently weak CIV lines 
associated with Lya lines of HI column density 10 14 B cm~ 2 
have been detected (Tytler 1995, Cowie et al. 1995) leading 
to [C/H] > -2.5 for those clouds. However no firm conclusion 
can be drawn from these detections. It shows that at least part 
of the Lya forest contains metal but the crucial observation 
would be to put strong limits on abundances at lower column 
density where the bulk of the Lya forest is to be found (Pe- 
titjean 1995a). Indeed the break in the TV(HI) column density 
distribution at 7V(HI) ~ 10 15 cm -2 has been interpreted as ev- 
idence for two populations of Lya systems (Petitjean et al. 
1993), the weakest (iV(HI) < 10 14 cm -2 ) and more numerous 
could be low metallicity intergalactic clouds and the strongest 
(7V(HI) > 10 15 cm~ 2 ) would be related to galactic haloes. In 
such a picture it is not surprising to find metals in the latter 
population. C IV and O VI doublets have been searched for in 
low-redshift Lya systems with no detection down to W r < 8 
and 10 mA respectively (Petitjean 1995b). 

Galaxies are strongly clustered and if the absorbing gas 
is associated with galaxies, absorption line systems are ex- 
pected to cluster on scales smaller than 1000 km s _1 . The issue 
has drawn much attention but is still unclear. At high redshift 
(z > 2), the distribution is practically uniform on scales larger 
than 300 km s _1 . On smaller scales, there has been some indi- 
cation for clustering of lines with W T > 0.36 X (Webb 1987, 
Ostriker et al. 1988, Crotts 1989, Barcons & Webb 1991). Since 
the bulk of the Lya forest has smaller equivalent width how- 
ever, high spectral resolution data should be used. Rauch et al. 
(1992) found no evidence for clustering using a sample of 295 
lines with 7V(HI) > 10 1375 cm -2 . Using similar data, Cher- 
nomordik (1995) claimed a detection of clustering on scales 
50-150 km s" 1 for lines with 7V(HI) > 10 14 cm" 2 . The fact 
that the properties may be column density dependent seems to 
be confirmed by analysis of recent new high resolution data 
(Cristiani et al. 1995). Great care has to be exercized however 
since the clustering signal can be diluted if the dimensions of the 
clouds are not negligeable compared to their mean separation 
(Bajtlik 1995). 

Sizes can be derived from detections of lines coincident in 
redshift in spectra of QSOs with small projected separations on 
the sky. The Lya forest appears strikingly similar in the two im- 
ages of the lensed quasar Q2345+007 separated by 7.3" (Foltz 
et al. 1984, Smette et al. 1992) indicative of large transversal 
dimensions. Only loose constraints are obtained with that ob- 
ject on the radius of the absorbers, R > 2/i _1 kpc. Stronger 



limits have been obtained from observation of brighter objects. 
HE1 104-1805A& B are separated by 3" and coincidences imply 
R > 50/t" 1 kpc for W T > 0.085 X (Smette et al. 1995). Bech- 
told et al. (1995) find 40/1- 1 < R < 280/j- 1 kpc at z ~ 1.8 
from Q1343+266A & B separated by 9.5". Dinshaw et al. 
(1995) derive 270 < R < 860/1^ kpc at z ~ 0.8 from FOS 
observation of 01 07-025 A & B separated by 1.44'. The latter 
dimensions are much larger than typical galactic sizes and rule 
out standard CDM minihalo models as well as models of freely 
expanding or pressure confined clouds. 

Association of Lya systems with galaxies is difficult to 
investigate directly, except at low redshift, where HST obser- 
vations are now available; galaxies are detected routinely using 
multi-object spectroscopy. Lanzetta et al. (1995) have found 
1 1 galaxies coincident in redshift with an absorption system 
of which nine are Lya-only systems and two have associated 
CIV absorptions. They report that (i) most of the z ~ 0.3 lu- 
minous galaxies are surrounded by gaseous envelopes of ra- 
dius ~160/i 10 1 kpc and (ii) the fraction of absorption clouds 
lying at a distance smaller than about \50h^ Q kpc from a lu- 
minous galaxy is at least 0.35±0.10 and likely 0.65±0.18. Le 
Brun et al. (1995) have detected 19 additional galaxies associ- 
ated with a Lya-only system and derived impact parameters as 
large as 500/ifgo kpc. However anti-correlation between Lya 
rest-frame equivalent width and the impact parameter is found 
only for the stronger lines. Moreover most of the studied Lya 
lines have W T > 0.30 X whereas the bulk of the Lya lines has 
W T < 0.3 X (Morris et al. 1991). The questions of whether 
weak absorptions are associated with galaxies or are completely 
unrelated to them (Morris et al. 1993) and what is the spatial 
distribution of the Lya gas are still open. 

In this respect, the search for Lya clouds in underdense 
regions and in voids of galaxies is of particular interest. At high 
redshift there is no clear evidence for regions devoid of lines 
(Carswell & Rees 1987, Duncan et al. 1989) although a few 
examples are known (Dobrzycki & Bechtold 1991, Rauch et al. 
1992). At low redshift, very recent HST observations of bright 
nearby AGNs by Stocke et al. (1995) have revealed one weak 
(W r = 36 mX) Lya line with no galaxy within 6 Mpc. 

In order to shed light on this issue, we have performed TV- 
body simulations of the kinematics of the dark matter in the 
universe, introducing a description of the warm photoionized 
gas. This gives a new framework in which observations can be 
interpreted. In a previous paper (Petitjean et al. 1995) we have 
shown that Lya gas traces the dark matter filaments. At low 
redshift, 25% of the lines should arise in extended envelopes 
of luminous galaxies of radius of the order of 1.5 Mpc, the 
remainder arising in gas independent of any galaxy and at a 
distance of 2 to 7 Mpc from the first neighbour galaxy. Here, 
we describe in detail our procedures and study the evolution 
with redshift of the number density of lines. Technical details 
about the simulations are given in Section 2. Results are given 
in Section 3 and conclusions are drawn in Section 4. 



2. Simulations 



We have improved and adapted the particle-mesh (PM) code 
described in detail by Kates et al. (1991) and Klypin & Kates 
(1991) to include the effects of photoionization. The PM code 
simulates the kinematics of the dark matter and the thermo- 
dynamics of baryons with periodic boundary conditions. The 
simulations reported here used 128 3 particles on a 256 3 grid. 
The choice of the length scales to be assigned to the simulation 
box and to the cells has to be considered carefully. At z = 0, 
structures are in the non-linear regime for scales smaller than 
10 Mpc (assuming h = 0.5). Unless otherwise stated, all results 
reported here were obtained using a simulation box of 12.8 
Mpc, corresponding to a comoving cell size of l c = 50 kpc. 
The physical cell size is of course l c /(l + z) 
The mass of one particle, M p , is equal to the mass of the matter 
distributed in a comoving volume of 8 cells with the mean 
background density today, p : 

M p = p (2l c ) 3 (1) 

If M p is smaller than the Jeans mass Mj, each particle can be 
considered as a non-collapsing, diffuse gas element, and the 
concentration of a sufficient large number of such elements de- 
scribes collapsing regions (filaments, merging processes etc..) 
Let us estimate the Jeans length for a gas with a baryonic mass 
density p gas and a temperature T, expressed in units of T 4 
(T = 10 4 T 4 ): 

where R is the gas constant, G the gravitational constant, /j, 
the molecular weight per hydrogen atom and fl h = < p gas > 
/p the average ratio of baryonic to dark-matter density. The 
requirement M p < Mj leads to the relation 

i 

2/ c < 1.201Mpc f — ) 2 (1 + z)~i (3) 

This relation is well satisfied up to high redshifts for our cell 
size. 

Our code averages the density field over the cell length. 
Let us regard one cell as an "elementary cloud." Let n\ = 
p/po denote the local density contrast. The maximum density 
contrast for an elementary cloud to be Jeans stable is then 
"i,max ~ (Aj,o/^ c ) 2 ~ 150 where Aj j0 is the Jeans length for 
p = p at z = 0. This density contrast corresponds to a HI 
column density iV Cjmax through the elementary cloud. For a 
nearly vanishing UV-background flux at z = one obtains 

\ 2 
j o 

-Wc.max « rii imax n H ,o'c = —r-no (4) 

where rin.o denotes the present hydrogen number density. 

For z > and an UV flux with an amplitude / in units of 
4.71-J21 with J 2 i = 10~ 21 erg s -1 sr -1 cm -2 it follows, 



« 7.7 10 15 /o _1 (^ir) "cm" 2 (6) 

Clouds with HI column densities smaller than 10 16 cm~ 2 at 
z = can thus be considered as Jeans stable. Clouds with 
higher column densities are unstable and should undergo frag- 
mentation. This happens only in the very central part of the dark 
matter halos. Higher resolution should be used to describe the 
structure of the densest regions. 

According to our procedure to be described below, the min- 
imum gas density for which a column density is computed is the 
mean baryonic density. This density corresponds to the smallest 
HI column density that the simulation can reproduce, A^ Cjmin . 
For z = 0, 7V Cjmin is smallerthan 10 11 cm~ 2 . For non-vanishing 
UV-background flux the z-dependence is given by 

N c ^(z) « °^ /o (7) 

It follows that for z > 3.3 and / = 0.5, the smallest HI 
column density to be found in a simulation with our resolution 
l c = 50 kpc is N(HT) ~ 10 13 cm" 2 . 

The positions, velocities, temperatures, and gas densities 
associated with the particles are followed starting from an ini- 
tial spectrum (details to be discussed below) and utilizing the 
Zel'dovich approximation as long as it remains valid. 

The baryonic material is approximated as following trajec- 
tories of the dark matter with constant Q b = 0.1. This approxi- 
mation is quite accurate in all regions prior to the first pancakes 
(or "shocks"). After shocking, the accuracy of the constant-f2 b 
approximation is of course degraded in regions where secondary 
shocking and/or very efficient cooling take place, but accuracy 
in these regions is not critical here. The conditions for maintain- 
ing constant fl h are well satisfied in Lya clouds: The potential 
wells of the dark matter structures must be deep enough to pre- 
vent the gas from escaping, but still shallow enough to keep the 
central matter density below levels permitting strong cooling 
and subsequent star formation to occur. 

As perturbations grow, the first objects start to collapse, pro- 
ducing shocks in the gas. Our code seeks shocks in the form of 
shell crossings. More precisely, a particle is labelled "shocked" 
if either of two criteria are met: (1) If the Jacobi determinant of 
the transformation from Lagrange (particle) to Euler (grid) co- 
ordinates becomes negative. (2) If the particle enters a shocked 
region, defined as entering a cell on a coarse grid of cell size 
2l c containing at least one shocked particle. In either case, we 
attempt to define a local velocity field U on the coarse grid, and 
if successful, we assign to the gas associated with the particle 
a temperature according to kT = ^ h wh(v - U) 2 /3, where v 
is the particle velocity, p,n the molecular weight and mn the 
mass of the hydrogen atom. (We take primordial abundances 
throughout.) It can happen in a small box that the tempera- 
ture assigned to the particle is below the minimum temperature 



(here 5000 K) to be considered "cold." In this case, we return 
the particle to the "unshocked" pool of particles, where it will 
most likely again be detected as shocked at some later timestep. 
The particle is also returned to the unshocked pool in those 
rare cases in which too few neighboring particles are present to 
determine U. The accuracy of this procedure in comparison to 
hydrodynamic simulations of pancake collapse is discussed in 
Kates etal. (1995). 

After successful assignment of temperature, we follow the 
thermal history by integrating the energy equation along particle 
trajectories in the form 



dT 
dt 



(7-1) 



T dun hm 1 



riH dt /xh krin 



(Mot — Tphot) 



(8) 



where 7 = 5/3 and rin is the number density of hydrogen atoms. 
Here, A to t is the total cooling rate (Klypin & Kates 1991, Kates 
et al. 1991) including radiative and dielectronic recombination, 
bremstrahlung, collisional excitation (Black 1981) and Comp- 
ton cooling (the latter is the most important process at z > 5). 

The heating rate due to photoionization by the ionizing 
background (i.e., from QSOs and galaxies) is denoted in Eq.(8) 
by r p hot- It is computed as the sum of the contributions by HI, 
Hel and Hell assuming ionization equilibrium 



£in(i) + £ e in(i)n e = n(i + l)n e a(i) 



(9) 



where i = (HI, Hel, Hell) and i + 1 = (HII, Hell, Helll), & 
and £ ei are the rates of the photo- and collisional ionization 
respectively, n e and n, the electronic and ionic number densi- 
ties, and a(i) denotes the recombination rate towards (i). The 
photo-ionization rate for each (i) is given by 



F v Oi{v) 



f 

6 ./, 



(10) 



where Iwi is the ionization energy (13.6 eV for HI). Opacity 
effects are taking into account writing 



F v = F 0v exp(-T„) 
with 



^rii(x)(j v (i) 



dx 



(11) 



(12) 



The photoionization cross-section of each of the species, a, de- 
creases with v faster than a v oc v~ 3 . Therefore we approximate 
the total optical depth t v by 77, (HI), T„(HeI) and T„(HeII) in 
the range v m < v < z/ He i, v HeI <v< v HeII and v HeII < v 
respectively. Therefore, 



&(x) =F eM-Ti(x))Gf 
where 

F 



G 
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Oi{y)dv 



(13) 



(14) 



The heating rate is written as 

F v (x)h(v - Vi 



f 



hv 



■Oi(y)dv 



F eM-n(x))Gfefn 



(15) 



where Gf and ef are respectively the ionization rate and heating 
paramter defined and given by Black (1981). 

To estimate the heating rate, we model each particle in 
accordance with the PM picture as an independent uniform 
cloud of volume (Z c ) 3 and density taken as the local density 
given by the code. Then we can assume Tj « erjrijZ c /2, which 
corresponds roughly to the optical depth near the center of this 
cloud. 

We solve equations (9) together with mass and charge con- 
servation and derive the heating and cooling rates assuming 
riHe = 0.1 riH and n e « nmi. The density ratios n(i+l)/n(i) 
are computed self-consistently, as described above. 
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Fig. 1. The UV flux generated by the code as function of 1 + z with 
/o = 0.57 (solid curve) and / = 0.78 (dashed curve). 



The shape of the ionizing spectrum is modelled as 
FvCxFoi/ -1 . F = /o47rJ2i is the ionizing flux at 13.6 eV 
and depends on redshift. It is assumed to be proportional to the 
rate, Am(T < 5000if ; z), at which baryonic material cools 
below 5000 K in the simulation: 



fo(z) = C cool Am(T 4 <0.5;z) 
+./;,(:+ A:/ 1 + ~ 



1 + z + Az 



(16) 



where C coo i is a factor of proportionality. This simulates an 
ionizing flux produced as a consequence of gas collapsing in the 



dense regions. In this case most of the photons just above 1 Ryd 
would have been produced by massive stars (see also Madau & 
Shull 1995). The redshift dependence is thus consistently given 
by the simulation; the only free parameter is the normalization 
(thus Ccooi) which is done assuming some value of / at z = 5. 
The flux as a function of 1 + z is given in Fig. 1 for two 
simulations with /o = 0.57 (solid curve) and / = 0.78 (dashed 
curve). In both cases, the flux remains nearly constant near its 
maximum value between z = 5 and z = 2.5. Below z = 2.5 the 
flux decreases approximately as a power law oc (1 + z)' 2 . The 
slope becomes steeper at z — > 0. 

The origin and intensity of the UV-background has been the 
subject of intensive investigations (e.g. Bechtold et al. 1987, 
Miralda-Escude & Ostriker 1990, Madau 1992). Although un- 
certainties are large and the interpretation questionable, mod- 
elling of the proximity effect gives a rough estimate of the inten- 
sity at high redshift (e.g. Bajtliketal. 1988,Luetal. 1991, Bech- 
told 1994) /o ~ 1 in the range 2 < z < 3. Somewhat smaller 
values (/ = 0.5), although depending on the shape of the ion- 
izing spectrum, have been derived from analysis of the statis- 
tical properties of metal line systems (Petitjean et al. 1992). 
At smaller redshift, only tentative estimates or upper limits 
are possible using the proximity effect (Kulkarny & Fall 1993, 
0.006 < /o < 0.04 at z ~ 0.5) or detection of Ha emission 
from nearby clouds (Songaila et al. 1989, Corbelli & Salpeter 
1993, Maloney 1993; 0.01 < / G < 0.1 at z ~ 0). Our results 
are consistent with what is expected for the time evolution. The 
value at low redshift is larger than the upper limit of Kulkarny 
& Fall (1993) by a factor of two but within the range from other 
determinations. 

3. Results 

3.1. Fluctuation spectrum 

The dynamics of large-scale structure formation is highly sensi- 
tive to the initial fluctuation spectrum. Confrontation of models 
with observations (see for example Gottlober et al. 1994 and 
Kates et al. 1995) strongly favor scenarios with broken scale 
invariance ("BSI"), mixed dark matter models, or low-mass 
models over "standard", COBE-normalized CDM. In partic- 
ular, CDM seems to imply an excess of virialized, massive 
(M > 10 13 M Q ) objects and too large velocity fields at small 
scales (as measured for example by relative pairwise velocities.) 
Nevertheless, one crucial property of COBE-normalized CDM 
convinced us to use it in the present simulations: the biasing 
factor of approximately unity. 

Most of the results reported here depend sensitively on an 
accurate representation of the statistics of the baryonic den- 
sity field. Now, as in any non-hydrodynamic simulation, the 
baryonic distribution is to be inferred from the dark-matter dis- 
tribution. Normally, if one is interested in large-scale structure 
and its observational consequences, a biasing factor may be in- 
troduced by some prescription, such as by searching for high 
peaks of the dark-matter distribution. This procedure roughly 
simulates the effects of baryonic condensations, and is thus rea- 



sonably well suited to a prediction of galaxy statistics. However, 
since it averages baryonic concentrations on a galactic scale, it 
does not appear to be suitable for an estimate of the density in 
Ly-a clouds. 

With a biasing factor of unity, we may safely assume a 
constant ratio of baryonic to dark-matter density without intro- 
ducing a systematic error in the baryonic density at the scales 
of interest. Within the limited spectral range resulting from the 
limited size of the computational box (see §2), the spectrum 
of most viable models differs from that of CDM by a nearly 
constant factor. The main features of CDM leading to an excess 
of massive structures will not be seen. Excessive small-scale 
power would admittedly be expected to influence our results 
somewhat, but certainly less than a systematic underestimate of 
the baryon density in clouds. Moreover, the amplitude of the 
power spectrum determines in general the epoch of the forma- 
tion of the first structures and thus the redshift dependence of 
the UV flux. Smaller amplitude implies later structure forma- 
tion. The UV peak is lower, but the distribution is shifted to 
later epochs, which can lead to a relatively higher background 
flux at small redshift. 

3.2. Probing the dark matter filaments 

From Fig. 1 one can see that the UV flux increases rapidly after 
z « 10. This rise corresponds to the redshift at which the first 
structures appear. 

Fig. 2 illustrates the evolution of the cloud distribution. 
With decreasing z, a trend toward concentration of the gaseous 
matter and a decrease in the number density is evident. The 
structures seen form the usual network of filaments connected 
by nodes where galaxies and groups of galaxies are to be found. 
The neutral hydrogen column density is computed through the 
slice, and contours of constant neutral hydrogen column density 
Af(HI) = 10 13 cm -2 are shown. It can be seen that the Lya gas 
traces the gravitational potential wells. 

3.3. Evolution of the line number density 

We now wish to compute the line number density evolution from 
our simulated, periodic universe. One copy of this universe 
consists of a box of comoving size L containing n 3 equally 
spaced cells and n 3 particles. For this computation, we imagine 
that a large number of light rays from ficticious, distant QSO's 
traverse our universe exactly once during a simulation timestep. 
(This condition fixes the timestep, of course). We therefore 
have a collection of many ray segments for each timestep. (A 
large number of rays is required for good statistics; in fact the 
number of rays exceeds the number of timesteps). Each ray is 
constructed by randomly choosing one ray segment per timestep 
from this collection. This procedure is intended to minimize the 
effects of periodicity and is in fact quite effective, as discussed 
in the following subsection. Rays intersect the x-y plane at cell 
vertices. 

Now suppose we have all the simulation data in the box at 
some time t , corresponding to a redshift z . We then perform 
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Fig. 2. Spacial distribution of Lya clouds in a 2 Mpc slice of a (25 
Mpc) 3 simulation box at redshifts z = (3.28, 1.6, 0.8) respectively. 
Gray contours show regions where N(H I) > 10 13 cm -2 through the 
box. 
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Fig. 3. ^-Dependence of the line number density dn/dz for lines with 
iV(HI) > 10 14 cm" 2 
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Fig. 4. «-Dependence of the line number density dn/dz for lines with 
N(m) > 10 14 cm -2 . Solid and dashed lines are for f = 0.57 and 
0.78 respectively. Observations are overplotted by open squares. 

one timestep of the simulation, obtaining the data at time t\ (red- 
shift zi). We may now construct the absorption line spectrum 
for redshifts between zq and z\ along the line of sight to our fic- 
ticious QSO: As mentioned above, the timestep At = t\ — to is 
taken as the time required for a light ray to traverse the comoving 
length L starting at zq. It turns out that At is always comfort- 
ably less than the maximum allowable timestep for numerical 
stability, i.e., particles move by less than a cell length within At. 
However, the time scale for temperature changes according to 
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Fig. 5. z-Dependence of of the line number density dn/dz for lines 
with iV(HI) > 10 13 cm~ 2 for two numerical resolutions, l c = 50 kpc 
(dotted line) and 25 kpc (solid line) 



Eq. (8) may be comparable to or shorter than At. We therefore 
integrate Eq. (8) along the ray with steps (Az) T « 10~ 3 A2:, 
where Az = Z\ — z . 

We define computational cells of 2 3 grid cells, i.e., n c /2 
per ray. We assign to each computational cell i along the line 
of sight a redshift z- x between zq and z\ (corresponding to a 
time t{) by interpolation as follows: The comoving distance 
traversed by the light ray between zq and Zi is simply i times 
the computational cell size 2l c , and z\ may then be obtained 
from 



2i« = f 

-no 



1 



1 



(17) 



where H denotes the Hubble constant and c the speed of light. 
Eq. (8) is then integrated for computational cell i up to ti. This 
procedure is iterated until the entire ray segment is covered by 
computational cells. With each computational cell, we associate 
the precise redshift Zi and all of the particle data (position, tem- 
perature, density) necessary for obtaining that computational 
cell's neutral hydrogen content and thus its contribution to the 
Lya line number density. Finally, we compile all the ray data 
for each ficticious QSO and store it for later analysis. 

We note for further reference that the resolution Sz of line 
contributions is clearly limited by the computational cell length: 
Szk2H I c (1 + z) 3 / 2 /c. 

The z-dependence of dn/dz is determined by averaging 
over our large sample of rays. Computation of dn/dz on any 
particular ray at some given redshift z = Z{ requires averaging 
over the contributions of all particle-clouds to the mean column 
density in computational cell i. We imagine that computational 
cell i contains many randomly spaced clouds with different 



temperatures and densities and that the distribution of sizes is 
restricted by the condition of pressure confinement. The av- 
eraging is thus to be weighted by the probability that the ray 
passes through a cloud. The mean column density N c (z) in one 
computational cell at z is 



N c (z) 



ENi(z) Wi 



(18) 



where the iV, (z) denote the column density of each particu- 
lar cloud in a computational cell at redshift position z. We recall 
that using a PM code each particle mass is distributed over a 
certain volume surrounding the particle situated in the center of 
this volume. In the case considered here this volume is equal to 
the volume of one grid cell (not computational cell). The over- 
lapping parts of these volumes define the mass contributions 
of nearby particles to the local density of the centered particle. 
Therefore, the "column density" of each cloud is formally equal 
to the density times the cell size. The w, in (18) are the proba- 
bilities that the ray crosses the cloud i. These probabilities are 
proportional to the cross section of the clouds, i.e. oc r\ where 
the r, are the radii of the clouds. The number density in a cloud 
is given as 



3m; 



Tli 



47TTO p r| 



(19) 



where mn denotes the mass of a hydrogen atom and nii = mo 
is the known mass of a cloud in the given simulation. Constant 
pressure leads using p = k^pTini = const, to the relation 



1 /3 

kT: ' , where k is a constant, and therefore to 



N c (z) 



£77 



(20) 



which is the mean column density contribution from the the 
computational cell under consideration on the simulated ray. 

We may define an effective cross section 7rr 2 ff of the com- 
putational cell as the average cross section of the particle clouds 
contained in the computational cell weighted by their column 
densities: 



'eff 



(21) 



The total volume (4tt/3)R 3 c of all particle clouds in the 



computational cell is £/, 



/t 3 £Tj. The corresponding 



cross section is 7ri? 2 c = k 2 7t(£ 7i) 2 / 3 . 

Therefore, the probability to find a line with column density 
N c (z) in the computational cell at redshift z on ray j is 



Wj(z) 



EjNjTf 1_ 

E, V, (£.7^)1 



(22) 



Finally, we estimate the line density dn/dz(z) for any given 
threshold column density N m by summing over all rays j those 



lines at redshift z whose column densities N c (z) exceed N H i, 
weighted by the probability Wj (z). 

Fig. 3 presents the results obtained for the z-dependence of 
the number density of clouds per unit redshift for clouds with 
-Whi > 10 14 cm~ 2 . The thick line presents the same results ob- 
tained by smoothing over redshift intervals of (A2;) smoo thing = 
0.2. For that line the following fit formulae can be given 



(IT) 

— = 16.8(1 + z) - 8 + 0.687(1 + z) 2 - 5 



(23) 



In Fig. 4 the line number density per unit redshift interval 
for column densities exceeding iVni = 10 14 cm -2 is plotted for 
the two different flux amplitudes used in Fig. 1 and compared 
with the data (square symbols). The latter are taken from com- 
pilations by Lu et al. (1991) and Petitjean et al. (1993) at high 
redshift and from HST observations by Bahcall et al. (1993). 
The computed number is also consistent with observation at 
z > 4 (Williger et al. 1994) although the latter are still scarce. 
The upper (solid) curve in Fig. 4 corresponds to the lower (solid) 
curve in Fig. 1. Obviously the upper curve, corresponding to 
/o = 0.57 yields a better fit to the observations. 

In Fig. 5, the line number density per unit redshift inter- 
val for column densities exceeding iVni = 10 13 cm~ 2 is plot- 
ted for the flux amplitude /o = .57. The thin solid line cor- 
responds to l c = 50 kpc. In the range from z = to 
z = 2.5, dn/dz behaves roughly like (1 + z) 0,7 . The ap- 
parent drop for z > 2.5 illustrates the redshift dependence of 
the resolution implied by Eq. (7). In order to extend the redshift 
range for iVni = 10 13 cm~ 2 , we carried out simulations with 
l c = 25 kpc. The result for corresponding / and the same 
TVhi limit is shown as the thick solid line in Fig. 5. To perform 
this simulation with higher resolution we had to decrease the 
box size. The scale comparable to that box size becomes non- 
linear already at redshift z « 1 . For this reason, the range of 
redshift shown is limited to z > z n \ > 1.5 . The observed num- 
ber of such lines is not well known yet due to the difficulties in 
observing weak lines at high enough resolution. However very 
recent observation with the Keck telescope show that at z ~ 3 
the number density of lines with N > 10 13 cm~ 2 is about 350 
per unit redshift (Songaila et al. 1 995). At low redshift the statis- 
tics of weak lines is very poor. From the Morris et al. ( 1 99 1 ) and 
Bahcall et al. (1991) data however it is possible to infer that the 
number per unit redshift of lines with W r > 25 mX is larger 
than 100 at z < 0. 15. Our computed numbers for / = 0.57 are 
a factor of two larger. However this can be fine-tuned using a 
somewhat larger / value. Given the limitations inherent to the 
simulation (see above), we consider that the discrepancy is not 
significant. It is very important to note that the bulk of the Lya 
forest has 7V(HI) < 10 14 cm~ 2 and that those lines evolve less 
rapidly than the stronger lines. The latter effect may be present 
in the data (Bechtold 1994, Bahcall et al. 1995). 

Fig 6 presents the column density distributions for a variety 
of cell lengths, redshifts, and fluxes: In general the curves repro- 
duce well the observed slope of (3 ss 1.7— 1.8for;z ss 3. In both 
curves representing z = 3 with low flux, a clear flattening is 
evident toward lower column densities. The flattening begins at 
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Fig. 6. Column density distribution dn/dN for assorted cell length, 
redshift, and flux values, solid line: l c = 50 kpc, z=0.5, /o = 0.57; 
dashed-multidotted line: l c = 50 kpc, z = 3 , /o = 0.57; dashed line: 
l c = 25 kpc, z = 3, fo = 0.57; dashed-dotted line: l c = 25 kpc, 
z = 3,/o = 1.1 
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7. Two-point correlation function £ c (ai) for clouds with 
= 10 14 cm" 2 (solid line) and for clouds with N c = 10 13 cm" 2 



(dashed line) at z=0. 
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Fig. 8. One dimensional correlation function for the absorption 
lines with Nhi > 10 14 cm~ 2 along the l.o.s. between z=3 and z=3 



higher TV in the l c = 50 kpc simulation than in the l c = 25 kpc 
simulation, consistent with the above remarks concerning res- 
olution. Higher flux (for l c = 25 kpc) decreases the flattening 
at low TV, in agreement with the prediction of Eq. (7). Further 
implications will be discussed in the Conclusions. 



3.4. Two point correlation 

Our computation is obviously insensitive to correlations on a 
spatial scale comparable to the box size; indeed, our ray segment 
procedure is designed to eliminate any spurious correlations 
resulting from the periodic box. Nevertheless, we may hope to 
detect correlations at some fraction of the box size: 

Fig. 7 shows the two-point correlation functions for 
clouds with TVhi = 10 14 cm~ 2 (solid line) and for clouds with 
TVhi = 10 13 cm -2 (dashed line) at z = 0. In both cases the cor- 
relation functions are bi-modal. A component with slope > - 1 is 
extended till separations of about 4 Mpc and indicates the con- 
tribution to filamentary and sheet-like structures. In the upper 
correlation function a spherical component is well expressed 
(slope < -2). The signal at small separations for the correlation 
function related to 7V H i = 10 14 cm -2 is much stronger than for 
the distribution of low density clouds. Both functions show a 
component with constant slope « -1.75 which is comparable 
with the slope of the two-point correlation function for the 
galaxy distribution. That might be understood as an indication 
of the association of those clouds with galaxies. 

In Fig. 8, the one-dimensional correlation function for 
the line distribution along a ray between z = 3 and z = 2 is 
shown for absorption lines with column densities exceeding 
10 14 cm~ 2 . The clear correlation out to a radius of about 4 
Mpc is in good accordance with the two-point correlation func- 



tion discussed above and with recent results of (Christiani et al. 
1995). The signal becomes weaker if we include the contribu- 
tion of low-density lines. 



4. Conclusions 

Our basic assumption was that the Ly-a clouds could trace the 
shallower potential wells given by the dark matter distribution. 
That concerns especially those dark matter structures as fila- 
ments and sheets. The quite good agreement of the obtained 
results with the observational data seems to indicate that the 
basis of our considerations is quite reasonable. 

On the other side this opens the exciting prospective to use 
Lya absorption line observations to map the universe in three 
dimensions and get information about the distribution of dark 
matter structures and their cosmological evolution. 

The weak dependence of the results with respect to the box 
size (provided it satisfies the initially posed restrictions) leads 
to the conclusion that the evolutionary effects in the number 
density of clouds are mainly due to the structure evolution on 
the smaller scales determining also the flux dependence on the 
redshift. 

Possible correlations of the clouds over scales larger than 
the box size cannot be investigated in the frame of our consider- 
ations due to the limited box size. However, the obtained suffi- 
ciently high signal in the one-dimensional correlation function 
for scales up to ss 5 Mpc has been proved with help of simula- 
tions with boxsize of about 25 Mpc. This result deserves further 
investigation and will be considered in a forthcoming paper. 

Incidentally, comparison of Ly-a data with sufficiently real- 
istic (i.e., hydrodynamic) simulations could be used in principle 
to "normalize" the power spectrum within the dynamical range. 
However, if one characterizes the "normalization" by an over- 
all amplitude, this amplitude will be overestimated in a small 
box due to truncation of long-wavelength modes. The relation 
between box size and time scales for structure formation (at 
constant amplitude) has been studied by Doroshkevich et al. 
(1995). 

Returning to the column density distribution, we saw that 
smaller cell length l c and higher flux f ss 1.0 reduce the 
flattening at low column densities. However, the data for 
iVffl > 10 14 cm~ 2 seem to favor a smaller flux / ss 0.6. In our 
simulations, the shielding of UV radiation in overdense regions 
is only taken into account within a cell. Our value of /o thus 
represents an underestimate of the flux that would be present in 
average or below average density regions (e.g., voids). A more 
accurate estimate of the behavior of the line distribution at low 
column densities would presumably need to take into account 
the distribution of small structures forming within "unshocked" 
regions, which were not included in our simulation. The contri- 
bution of such structures within voids could be estimated using a 
"mini-halo" model in combination with Press-Schechter theory 
(Mo, Miralda-Escude, & Rees, 1993). 
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